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Abstract 


A stochastic background of gravitational waves is expected to arise from a superposition of many 
incoherent sources of gravitational waves, of either cosmological or astrophysical origin. This background 
is a target for the current generation of ground-based detectors. In this article we present the first joint 
search for a stochastic background using data from the LIGO and Virgo interferometers. In a frequency 
band of 600-1000 Hz, we obtained a 95% upper limit on the amplitude of Hqw (/) = ^3 (//900Hz) 3 , of 
0,3 < 0.33, assuming a value of the Hubble parameter of h \ oo = 0.72. These new limits are a factor of 
seven better than the previous best in this frequency band. 
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I. INTRODUCTION 


A major science goal of current and future generations of gravitational- wave detectors is the de- 
tection of a stochastic gravitational wave background (SGWB) - a superposition of unresolvable 
gravitational-wave signals of astrophysical and/or cosmological origin. An astrophysical back- 
ground is expected to be comprised of signals originating from astrophysical objects, for example 
binary neutron stars [1], spinning neutron stars [2], magnetars [3] or core-collapse supemovae [4]. 
A cosmological background is expected to be generated by various physical processes in the early 
universe [5] and, as gravitational waves are so weakly interacting, to be essentially unattenuated 
since then. We expect that gravitational waves would decouple much earlier than other radiation, 
so a cosmological background would carry the earliest information accessible about the very early 
universe [6]. There are various production mechanisms from which we might expect cosmological 
gravitational waves including cosmic strings [7], amplification of vacuum fluctuations following 
inflation [8, 9], pre-Big-Bang models [10, 1 1], or the electroweak phase transition [12]. 

Whatever the production mechanism of a SGWB, the signal is usually described in terms of the 
dimensionless quantity, 

sw/) = 4^, CD 

pc dj 

where dpow is the energy density of gravitational radiation contained in the frequency range / to 
f + df and p c is the critical energy density of the universe [13]. As a SGWB signal is expected 
to be much smaller than current detector noise, and because we assume both the detector noise 
and the signal to be Gaussian random variables, it is not feasible to distinguish the two in a single 
interferometer. We must therefore search for the SGWB using two or more interferometers. The 
optimal method is to cross-correlate the strain data from a pair, or several pairs of detectors [5, 13]. 
In recent years, several interferometric gravitational wave detectors have been in operation in the 
USA and Europe. At the time that the data analysed in this paper were taken, five interferometers 
were in operation. Two LIGO interferometers were located at the same site in Hanford, WA, one 
with 4km arms and one with 2km arms (referred to as HI and H2 respectively). In addition, one 
LIGO 4km interferometer, LI, was located in Livingston, LA [14]. The Virgo interferometer, VI, 
with 3km arms was located near Pisa, Italy [15] and GEO600, with 600m arms, was located near 
Hannover, Germany [16]. LIGO carried out its fifth science run, along with GEO600, between 
5th November 2005 and 30th September 2007. They were joined from 18th May 2007 by Virgo, 
carrying out its first science run. In this paper we present a joint analysis of the data taken by 
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the LIGO and Virgo detectors during these periods, in the frequency range 600-1000 Hz. This 
is the first search for a SGWB using data from both LIGO and Virgo interferometers, and the 
first using multiple baselines. Previous searches using the LIGO interferometers used just one 
baseline. The most sensitive direct limit obtained so far used the three LIGO interferometers, but 
as the two Hanford interferometers were colocated this involved just one baseline [17]. The most 
recent upper limit in frequency band studied in this paper was obtained using data from the LIGO 
Livingston interferometer and the ALLEGRO bar detector, which were colocated for the duration 
of the analysis [18]. The addition of Virgo to the LIGO interferometers adds two further baselines, 
for which the frequency dependence of the sensitivity varies differently. The frequency range used 
in this paper was chosen because the addition of Virgo data was expected to most improve the 
sensitivity at these high frequencies. This is due in part to the relative orientation and separation 
of the LIGO and Virgo interferometers, and in part to the fact that the Virgo sensitivity is closest 
to the LIGO sensitivity at these frequencies. The GEO600 interferometer was not included in this 
analysis as the strain sensitivity at these frequencies was insufficient to significantly improve the 
sensitivity of the search. 

The structure of this paper is as follows. In Section II we describe the method used to anal- 
yse the data. In Section III we present the results of the analysis of data from the LIGO and 
Virgo interferometers. We describe validation of the results using software injections in Section 
IV. In Section V we compare our results to those of previous experiments and in Section VI we 
summarise our conclusions. 


II. ANALYSIS METHOD 


The output of an interferometer is assumed to be the sum of instrumental noise and a stochastic 
background signal, 

s(t) = n(t ) + h[t). (2) 

The gravitational wave signal has a power spectrum, Sgw(/)> which is related to Hgw(/) by [19] 

3Hq n GW (/) 


Scw(f) = 


Our signal model is a power law spectrum, 


IOtj-2 p 


f 


(3) 


^GW (/) = ( ~T ) , 

JR , 


(4) 
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where a is the spectral index, and f R a reference frequency, such that Q n = Qgw (//,')■ For this 
analysis we create a filter using a model which corresponds to a white strain amplitude spectrum 
and choose a reference frequency of 900 Hz, such that 


^Gw(/) 




/ 

900Hz 


(5) 


We choose this spectrum as it is expected that some astrophysical backgrounds will have a rising 
9 gw (/) spectrum in the frequency band we are investigating [2-4]. In fact, different models 
predict different values of the spectral index a in our frequency band, so we quote upper limits for 
several values. 

For a pair of detectors, with interferometers labelled by i and j, we calculate the cross- 
correlation statistic in the frequency domain 


Y = 


/*+ OO 

' —oo 
/* + OO 


df Y (/) 


r»+oo 


df 


df StU - f) ~s*(f) ~ Sj (f) Qij(f ) , 


(6) 


where S;(/) and sf f ) are the Fourier transforms of the strain time-series of two interferometers, 
Qij{f) is a filter function and dj is a finite-time approximation to the Dirac delta function, [13] 


Sr(f) ■= 


r dfe-^< = Sin(,r / r) 

l-T/2 71 J 


(7) 


We assume the detector noise is Gaussian, stationary, uncorrelated between the two interferometers 
and much larger than the signal. Under these assumptions, the variance of the estimator Y is 


a 


2 

Y 



df Oy(f) 


ry t- + OO 

« 2 J o dfPi(f)Pff) I Qij(f) I 2 , ( 8 ) 

where Pff) is the one-sided power spectral density of interferometer i and T is the integration 
time. By maximizing the expected signal-to-noise ratio (SNR) for a chosen model of Hgw (/), we 
find the optimal filter function, 


QaU)= V 


7 «(/)/ 


Q — 3 


(9) 


m(mv) 1 

where 7 jj (/) is the overlap reduction function (ORF) of the two interferometers and AT is a normal- 
isation factor. We choose the normalisation such that the cross-correlation statistic is an estimator 
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of O a , with expectation value (Y) = Q n . It follows that the normalisation is 



( 10 ) 


Using this filter function and normalisation gives an optimal SNR of [13] 



(ID 


The ORF encodes the separation and orientations of the detectors and is defined as [13, 20] 


is the response of the ith detector to the A = +, x polarisation, where e^ b are the transverse 
traceless polarisation tensors. The geometry of each interferometer is described by a response 
tensor, 


which is constructed from the two unit vectors that point along the arms of the interferometer, x 
and y [20, 21]. At zero frequency, the ORF is determined solely by the relative orientations of 
the two interferometers. The LIGO interferometers are oriented in such a way as to maximize the 
amplitude of the ORF at low frequency, while the relative orientations of the LIGO- Virgo pairs 
are poor. Thus at low frequency the amplitude of the ORF between the Hanford and Livingston 
interferometers, 7 HL(f), is larger than that of the overlap between Virgo and any of the LIGO 
interferometers, 7 Hv(f) or 7 Lv(f) (note that the ‘HL’ and ‘HV’ overlap reduction functions hold 
for both HI and H2 as they are colocated). However, at high frequency the ORF behaves as a sine 
function of the frequency multiplied by the light travel time between the interferometers. As the 
LIGO interferometers are closer to each other than to Virgo, their ORF 7 mif) oscillates less, but 
decays more rapidly with frequency than the the ORFs of the LIGO-Virgo pairs. Fig. 1 shows the 
ORFs between the LIGO Hanford, LIGO Livingston and Virgo sites. 

We define the “sensitivity integrand”, 1 (f), by inserting Eq. 9 and Eq. 10 into Eq. 8 , giving 



( 12 ) 


where His a unit vector specifying a direction on the two-sphere, Ax = x % — x 3 is the separation 
of the two interferometers and 



(13) 



(14) 



( 15 ) 
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f (Hz) 



Figure 1 : Plot of the overlap reduction function (ORF) for the pairs of sites used in this analysis. The dashed 
curve is the ORF for the two LIGO sites (F1L), the solid curve is for the Flanford-Virgo sites (HV) and the 
dashed-dotted curve is for Livingston- Virgo (LV). We see that the LIGO orientations have been optimized 
for low frequency searches, around 10-100 Hz. However, this ORF falls off rapidly with frequency, such 
that at frequencies over ~ 500 Hz, the amplitude of the ORF of the HL pair is smaller than that of the Virgo 
pairs. The LV and HV overlap reduction functions oscillate more with frequency, but fall off more slowly, 
due to the larger light-travel time between the USA and Europe. 

This demonstrates the contribution to the inverse of the variance at each frequency. The sensitivity 
of each pair is dependent on the noise power spectra of the two interferometers, as well as the 
observing geometry, described by 7 For interferometers operating at design sensitivity, this 
means that for frequencies above ~ 200 Hz the LIGO-Virgo pairs make the dominant contribu- 
tion to the sensitivity [22] . During its first science run Virgo was closest to design sensitivity at 
frequencies above several hundred Hz, which informed our decision to use the 600-1000 Hz band. 

The procedure by which we analysed the data is as follows. For each pair of interferometers, 
labelled by I, the coincident data were divided into segments, labelled by J, of length T = 60s. 
The data from each segment are Hann windowed in order to minimize spectral leakage. In order 
not to reduce the effective observation time, the segments are therefore overlapped by 50%. For 
each segment, the data from both interferometers were Fourier transformed then coarse-grained 
to a resolution of 0.25 Hz. The data from the adjacent segments were then used to calculate 
power spectral densities (PSDs) with Welch’s method. The Fourier transformed data and the 
PSDs were used to calculate the estimator on f2 3 , V}./, and its standard deviation, a rj . For each 
pair, the results from all segments were optimally combined by performing a weighted average 
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(with weights 1 /erf j), taking into account the correlations that were introduced by the overlapping 
segments [23]. The weighted average for each pair, Y r , has an associated standard deviation 
ay, also calculated by combining the standard deviations from each segment (note that ay is the 
equivalent of cry (from Eq. 8) for each pair, /, but we have dropped the Y subscript to simplify the 
notation). 


A. Data Quality 


Data quality cuts were made to eliminate data that was too noisy or non- stationary, or that 
had correlated noise between detectors. Time segments that were known to contain large noise 
transients in one interferometer were removed from the analysis. We also excluded times when 
the digitizers were saturated, times with particularly high noise and times when the calibration 
was unreliable. This also involved excluding the last thirty seconds before the loss of lock in the 
interferometers, as they are known to have an increase in noise in this period. Additionally, we 
ensured that the data were approximately stationary over a period of three minutes, as the PSD 
estimates, Pi(f), used in calculating the optimal filter and standard deviation in each segment are 
obtained from data in the immediately adjacent segments. This was achieved by calculating a 
measure of stationarity, 


Actij = 


\<TIJ 


a 


ij I 


<Jij 


(16) 


for each segment, where oyj was calculated (following Eq. 8) using the PSDs estimated from the 
adjacent segments, and a' LJ was calculated using the PSDs estimated using data from the segment 
itself. To ensure stationarity, we set a threshold value, (, and all segments with values of A 07 j > ( 
are discarded. The threshold was tuned by analysing the data with unphysical time offsets between 
the interferometers; a value of ( = 0.1 as this ensures that the remaining data are Gaussian. 

In order to exclude correlations between the instruments caused by environmental factors we 
excluded certain frequencies from our analysis. The frequency bins to be removed were identified 
in two ways. Some correlations between the interferometers were known to exist a priori, e.g. there 
are correlations at multiples of 60 Hz between the interferometers located in the USA due to the 
frequency of the power supply [19]. These were removed from the analysis, but in order to ensure 
that all coherent bins were identified, we also calculated the coherence, 


r(/) = 


|gi(/)g 2 (/)| 2 

Pitmif) ’ 


(17) 
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which is the ratio of the cross-spectrum to the product of the two power spectral densities, averaged 
over the whole run. This value was calculated first at a resolution of 0. 1 Hz, then at 1 mHz to 
investigate in more detail the frequency distribution of the coherence. Several frequencies showed 
excess coherence; some had been identified a priori but two had not, so these were also removed 
from the analysis. The calculations of the power spectra and the cross-correlation were carried out 
at a resolution of 0.25 Hz, so we removed the corresponding 0.25 Hz bin from our analysis. Excess 
coherence was defined as coherence exceeding a threshold of T(/) = 5 x 1 0 “ 3 . This threshold 
was also chosen after analysing the data with unphysical time offsets. The excluded bins for each 
interferometer can be seen in Table I. 


IFO 

Notched frequencies (Hz) 

HI 

786.25 

Harmonic of calibration line 


961 

Timing diagnostic line 

H2 

640 

Excess noise 


814.5 

Harmonic of calibration line 


961 

Timing diagnostic line 

LI 

793.5 

Harmonic of calibration line 


961 

Timing diagnostic line 

VI 

706 



710 

Harmonics of 


714 

calibration lines 


718 



Table I: Table of the the frequency bins excluded from the analysis for each interferometer. The bins at 640 
Hz and 961 Hz were identified using coherence tests, while the others were excluded a priori. Also excluded 
were harmonics of the power line frequency at multiples of 60 Hz for the LIGO detectors and multiples of 
50 Hz for Virgo. Each excluded bin is centred at the frequency listed above and has a width of 0.25 Hz. 


B. Timing accuracy 

In order to be sure that the cross-correlation is a measure of the gravitational-wave signal 
present in both detectors in a pair, we must be sure that the data collected in both detectors are truly 
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coincident. Calibration studies were carried out to determine the timing offset, if any, between the 
detectors and to estimate the error on this offset. These studies are described in more detail in 
reference [24], but we summarize them here. 

The output of each interferometer is recorded at a rate of 16384 Hz. Each data point has 
an associated time-stamp and we need to ensure that data taken with identical time stamps are 
indeed coincident measurements of the strain, to within the calibration errors of the instruments. 
No offset between the instruments was identified, but several possible sources of timing error 
were investigated. First, approximations in our models of the interferometers can introduce phase 
errors. For the measurement of strain, we model the interferometers using the long- wavelength 
approximation (i.e. we assume that the wavelengths of the gravitational waves that we measure are 
much longer than the arm-lengths of the interferometers). We also make an approximation in the 
transfer function of the Fabry-Perot cavity; the exact function has several poles or singularities, 
but we use an approximation which includes only the lowest frequency pole [25]. The errors that 
these two approximations introduce largely cancel, with a residual error of ~ 2/is or ~ 1° at 1kHz 
[24], 

Secondly, there is some propagation time between strain manifesting in the detectors and the 
detector output being recorded in a frame file. This is well understood for all detectors and is 
accounted for (to within calibration errors) when the detector outputs are converted to strain. The 
time-stamp associated with each data point is therefore taken to be the GPS time at which the 
differential arm length occurs, to within calibration errors [24]. 

Thirdly, the GPS time recorded at each site has some uncertainty. The timing precision of the 
GPS system is ~ 30 ns, which corresponds with the stated location accuracy of ~ 10 m. Each site 
necessarily uses its own GPS receiver, so the relative accuracy of these receivers has been checked, 
by taking a Virgo GPS receiver to a FIGO site and comparing the outputs. The relative accuracy 
was found to be better than 1/is. The receivers have also been checked against Network Time 
Protocol (NTP) and were found to have no offset [24] . The total error in GPS timing is far smaller 
than the instrumental phase calibration errors in the 600-1000 Hz frequency band (see Table II). 

These investigations concluded that the timing offset between the between the instruments is 
zero for all pairs, with errors on these values that are smaller than the error in the phase calibration 
of each instrument. The phase calibration errors of the instruments are negligible in this analysis 
as their inclusion would produce a smaller than 1% change in the results at this sensitivity, and 
therefore the relative timing error is negligible. 
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C. Combination of multiple pairs 


We performed an analysis of all of the available data from LIGO’s fifth science run and Virgo’s 
first science run. However, we excluded the H1-H2 pair as two instruments were built inside 
the same vacuum system, and so may have significant amounts of correlated noise. There is an 
ongoing investigation into identifying and removing these correlations [26], and for the present 
analysis, we consider only the five remaining pairs. As described above, the output of each pair 
yields an estimator, Yj, with a standard deviation, 07 , where I = 1 ... 5 labels the detector pair. 

Using the estimators Yj and their associated error bars, 07 , we construct a Bayesian posterior 
probability density function (PDF) on 0 3 . Bayes theorem says that the posterior PDF of a set of 
unknown parameters, 6, given a set of data, V, is given by 


p(0\V) 


Pimm 

p{V) 


(18) 


where p(6) is the prior PDF on the unknown parameters - representing the state of knowledge 
before the experiment -p(T>\9) is the likelihood function - representing the probability distribution 
of the data given particular values of the unknown parameters - and p(V) is a normalisation factor. 
In this case, the unknown parameters, 9, are the value of and the amplitude calibration factors 
of the instruments, which will be discussed below. The data set, V, is the set of five estimators, 
{Yj}, we obtain from the five pairs of instruments. 

In forming this posterior, we must consider the errors in the calibration of the strain data ob- 
tained by the interferometers. In the data from one interferometer, labelled by i, there may be an 
error on the calibration of both the amplitude and the phase, such that the value we measure is 


*(/) = 


_ A 


S'iif), 


(19) 


where s*(/) is the “true” value that would be measured if the interferometer were perfectly cali- 
brated. The phase calibration errors given in Table II are negligible, and the studies described in 
Section II B have shown that there is no significant relative timing error between the interferom- 
eters, so we can simply assume that 0, = 0. However, the amplitude calibration errors are not 
negligible, and the calibration factors take the values A* = 0 ± e\ 7 , where e \ 7 are the fractional 
amplitude calibration errors of the instruments, which are quoted in Table II. 

The calibration factors combine such that the estimator for a pair / is 


Yj = e Al ’ 1+Al ’ 2 Yj , 


( 20 ) 
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Instrument 

Amplitude error (%) 

Phase error (deg) 

HI 

10.2 

4.3 

H2 

10.3 

3.4 

LI 

13.4 

2.3 

VI 

6.0 

4.0 


Table II: Table of values of the errors in the calibration of amplitude and phase for each of the LIGO [27] 
and Virgo [28] instruments used in this analysis. The errors are valid over the whole 600-1000 Hz band. 


where Yf is the “true” value that would be measured with perfectly calibrated instruments and A/ , 
and A/ 2 are the calibration factors of the two instruments in pair I. The likelihood function for a 
single estimator is is given by 


p(y7|G 3 ,(T/, A/) 

where we have used A T = A I}1 + A/ 2 
over all pairs of Equation 21 


/ 


: exp 


(y t - e A/ G 3 ) 


! \ 


V 


2 erf 


( 21 ) 


cr/ \/2vr 

. The joint likelihood function on all the data is the product 


/ 


^-pairs 

p({Y r }|0 3 ,{(T/},{A / }) = Yl crj, A/). 

i=i 


( 22 ) 


In order to form a posterior PDF, we define priors on the calibration factors of the individual in- 
terferometers, {A,}. The calibration factors are assumed to be Gaussian distributed, with variance 
given by the square of the calibration errors quoted in Table II, such that 


niFO 1 / a 2 \ 

PdMIKd) = II — -f , (23) 

i— i GV,iV 27T y ,i J 

where tiifo is the number of interferometers we are using, in this case four. The prior on G 3 is a 
top hat function 


P(tt S) 


q for 0 ^ <C f^max 

‘max 


(24) 


I 0 otherwise 

We choose a flat prior on 0 3 because, although there has been an analysis in this band previously, 
it did not include data from the whole of the frequency band and an uninformative flat prior is 
conservative. We chose Q max = 10, which is two orders of magnitude greater than the estimators 
and their standard deviations, such that the prior is essentially unconstrained. 
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We combine the prior and likelihood functions to give a posterior PDF 


p(D 3 ,{Ai}|{yf},{a/},{6 A ,i}) = p(D 3 )p({Ai}|{e A ,i}) . . . (25) 

x P ({y>}|n 3 ,W,{A,}). 

We marginalize this posterior analytically over all A, to give us a posterior on Q 3 alone, 


p(D 3 |{Y/},{oy},{e Ai i}) = / dA 1 dA 2 . . . / dA niFO p(Q 3 , {Ai}|{Y}}, {oy}, {e Aji ]#6) 


Using this posterior PDF we calculate a 95% probability interval, (f2i ower , f2 uppe r) on Q 3 . We 
calculate the values of fli ower and O upper by finding the minimum-width interval that satisfies 

P ^UDDer 


^lo 


p(fi 3 |{y>}, {tr/}, {e Aii })dU 3 = 0.95. 


(27) 


If we find that 0 lower is equal to zero, then we have a null result, and we can simply quote the 
upper limit, O upper . 

The optimal estimator, Y, is given by the combination of Yj, a j and A j that maximises the 
likelihood, such that 

Ei^aJ 2 


Y = 




It has a variance, cr, given by 


a ~ 2 = 


E 


e 2A ^7 2 . 


(28) 


(29) 


Under the assumption that the calibration factors A* are all equal to zero, then the optimal way to 
combine the results from each pair is to perform a weighted average with weights 1/crf (equiva- 
lently to combining results from multiple, uncorrelated, time segments) [22] 


Y = 


a~ 2 = 


E^J 2 


E 


-2 


CTr 


(30) 

(31) 


A combined sensitivity integrand can also be found by summing the integrands from each pair: 

[ 22 ] 

(32) 
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III. RESULTS 


We applied the analysis described in Section II to all of the available data from the LIGO and 
Virgo interferometers between November 2005 and September 2007 1 and obtained estimators of 
0 3 from each of five pairs, which are listed in Table III along with their standard deviations. We 
also create the combined estimators and their standard deviations, using Equations 30 and 31, for 
the full network, and for the network including only the LIGO interferometers. We see that the 
addition of Virgo to the network reduces the size of the standard deviation by 23%. 


Network 

Estimator Yj 

H1L1 

0.11 ±0.15 

HIV 1 

0.57±0.22 

H2L1 

-0.14±0.26 

H2V 1 

-0.52±0.41 

L1V 1 

0.18±0.20 

LIGO 

0.05±0.13 

all 

0.15±0.10 


Table III: Table of values of Yj, the estimator of 0 3 , obtained by analysing the data taken during LIGO’s fifth 
science run and Virgo’s first science run, over a frequency band of 600-1000 Hz, along with the standard 
deviation, 07 , of each result. 

Using the posterior PDF defined in Equation 26 and the calibration errors in Table II we found 
a 95% upper limit of 0 3 < 0.33, assuming the Hubble constant to be /ii 00 = 0.72 [29], while 
using only the LIGO instruments obtained an upper limit of 0 3 < 0.31. Both of the lower limits 
were zero. The posterior PDFs obtained by the search are shown in Figure 2, while the sensitivity 
integrands, which show the contribution to the sensitivity of the search from each frequency bin, 
are shown in Figure 3. The upper limit corresponds to a strain sensitivity of 8.5 x 10- 24 Hz~ 1/2 
using just the LIGO interferometers, or 8.7 x 10 24 Hz 1/2 using both LIGO and Virgo. The 
LIGO-only upper limit is, in fact, lower than the upper limit using the whole data set, even though 

1 We initially analyzed only data from times after Virgo had begun taking data (May-September 2007). This pre- 
liminary analysis resulted in a marginal signal with a false-alarm probability of p=2%. To follow up, we extended 
the analysis to include all available LIGO data, yielding the results shown here, which are consistent with the null 
hypothesis. 
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the sensitivity of the combined LIGO-Virgo analysis is better. This is not surprising because 
the addition of Virgo also increases the value of the estimator. The estimator will usually lie 
somewhere between 0 and 2 a - in this case, the LIGO-only estimator was in the lower part of that 
range while the LIGO-Virgo estimator was not, but the two results are entirely consistent with each 
other. When we add Virgo, the likelihood excludes more of the parameter space below Q : > = 0, 
but this is a region we already exclude by setting the priors. Monte-Carlo simulations show that, 
in the absence of a signal, the probability of the combined LIGO-Virgo upper limit being at least 
this much larger than the LIGO-only upper limit is 4.3%. This probability is not so small as to 
indicate a non-null result and we therefore conclude that the LIGO-Virgo upper limit is larger due 
to statistical fluctuations. 

We also used the same data to calculate the 95% probability intervals for gravitational wave 
spectra with spectral indices ranging over —4 < a < 4, which correspond with different models 
of possible backgrounds in our frequency band. For example, a background of magnetar signals 
would be expected to have a spectral index of a = 4 [3]. Figure 4 shows the values of these upper 
limits. Note that they were all calculated using a reference frequency of 900 Hz, and Hubble 
parameter h wo = 0.72. 



Figure 2: Posterior PDFs on Os.The dashed line shows the posterior PDF obtained using just the LIGO 
detectors, the solid line shows the PDF obtained using LIGO and Virgo detectors. The filled areas show the 
95% probability intervals. 
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Figure 3: Sensitivity integrands for the LIGO only result (dashed) and for the full LIGO- Virgo result (solid). 
We can see that the sensitivity is increased across the band by the addition of the Virgo interferometer to the 
search. The vertical lines correspond to frequency bins removed from the search. 



Figure 4: 95% probability intervals on Q a , calculated using different values of a. These upper limits were 
all calculated using the same data, with a band width of 600-1000 Hz and a reference frequency of 900 Hz. 
The dashed line shows the upper limit calculated using the LIGO interferometers only, while the solid line 
shows the upper limits calculated using all of the available data. The lower limits were all zero. 

IV. VALIDATION OF RESULTS 

In order to test our analysis pipeline, we created simulated signals and used software to add 
them to the data that had been taken during the first week of Virgo’s first science run (this week was 
then excluded from the full analysis). We generated frame files containing a simulated isotropic 
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Figure 5: Plot of the recovered values of Q :j for five software injections. The error bars show the 95% 
probability intervals. The quietest injection had a lower limit equal to zero. Note that all analyses excluded 
H2. Each injection used the same data, with the simulated signal scaled to different amplitudes. 

stochastic background, with 0 GW (/) oc / 3 . We were then able to scale this signal to several values 
of fl 3 and add it to the data taken from the instruments. We did not include H2 in this analysis, but 
used only HI, LI and VI. Table IV shows the injected values of fl 3 and the recovered values and 
associated standard deviations, along with the SNR of the signal in the HI VI pair. The recovered 
95% probability intervals of the injections can be seen in Figure 5. The intervals all contain the 
injected value of fi 3 . 

It should be noted that, in order to have detectable signals in this short amount of data, the larger 
injections are no longer in the small signal limit. We usually make two assumptions based on this 
limit. The first is the approximation in Eq. 8, which only holds if the signal is much smaller than 
the noise, as we are ignoring terms that are first and second order in Q. G \v(/) [13]. The second as- 
sumption enters into the calculation of the noise PSDs, Pi(f)- We calculate these directly from the 
data, as in the small-signal limit we can assume that (|sj(/)| 2 ) ~ (|nj(/)| 2 ).The first assumption 
causes an over-estimation of the standard deviation, while the second causes our “optimal” filter to 
no longer be quite optimal. If we ignore these assumptions, we will underestimate the theoretical 
error bar, cry , and the width of the posterior PDFs. However, we still find 95% probability intervals 
that are consistent with the injected signals. 
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Injected H 3 

Estimator Y 

95% probability interval 

SNR in HIV 1 

2.0 

1.8T1.3 

(0.0, 4.1) 

1.3 

9.7 

9.1T1.5 

(5.7, 12.8) 

6.3 

20.2 

19.3T1.8 

(14.2,24.8) 

13.3 

95.1 

91.1T3.7 

(72.3, 110.6) 

62.3 

203.1 

194.1T6.2 

(154.9,234.3) 

133.1 


Table IV: Table of values of Q 3 for software injections, along with the recovered values, the 95% proba- 
bility interval and the expected SNR of each injection in the HI VI pair. Note that the standard deviations 
presented in this table are underestimated, as the injections are not in the small signal limit, however we 
still recover the signals within the 95% probability intervals. 

V. COMPARISON WITH OTHER RESULTS 

The previous most sensitive direct upper limit in this frequency band was Ogw (/) < 1-02, 
obtained by the joint analysis of data from the LIGO Livingston interferometer and the ALLEGRO 
bar detector over a frequency band of 850 Hz < / < 950 Hz [18]. This result was obtained using 
a constant fi G w (/) — H 0 , so should be compared with our upper limit for a = 0. As can be seen 
in Figure 4, our 95% upper limit for a = 0 is < 0.16 using all the available data, or < 0.15 
using just the LIGO interferometers, therefore our result has improved on the sensitivity of the 
LIGO-ALLEGRO result by a factor of ~ 7. The comparative strain sensitivity of the upper limits 
of the current search and the LIGO-ALLEGRO search can be seen in Figure 6. 

The previous most sensitive direct limit at any frequency was the analysis of data from the three 
LIGO detectors in the fifth science run [17]. The analysis was carried out using the same data as 
the analysis presented in this paper, but was restricted to the frequency band 40 Hz < / < 500 Hz. 
This included the most sensitive frequency band of the three detectors. The 95% upper limit on 
ff 0 in this band was given as 6.9 x 10 -6 , which is a factor of 2 x 10 4 times smaller than our 
upper limit. They also found an upper limit on Q :J of 7.1 x 10 -6 . In order to compare that to our 
upper limit on fl 3 , we must extend the spectrum to the frequency band analysed in this paper. The 
40 Hz < / < 500 Hz upper limit would correspond to an upper limit at 900 Hz of fi 3 < 0.0052, 
which is a factor of ~ 60 smaller than the upper limit presented in this paper. The search at lower 
frequencies is significantly more sensitive and we would expect that in the advanced detector era 
the combined analysis of LIGO and Virgo detectors at low frequencies will improve even further 
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Figure 6: Comparison of the strain sensitivity of two searches for an isotropic stochastic background of 
gravitational waves. The two solid grey lines show strain sensitivity of the Hanford 4km interferometer 
(dark grey) and the Virgo interferometer (light grey), these spectra were obtained by averaging over the 
data analysed in this paper. The dot-dashed line shows the main result of this paper, the search for a 
SGWB with Hgw(/) °c / 3 , which is white in strain amplitude, and corresponds to an upper limit of 
Q 3 < 0.33. The dashed line shows the result of the same search, but for constant Hgw(/)> and corresponds 
to an upper limit of Ho < 0.16. The solid black line shows the strain sensitivity of the LIGO-ALLEGRO 
search, which corresponds to an upper limit of Ho < 1.02 and was calculated over a frequency range of 
850 Hz < / < 950 Hz [18]. The two dotted lines show the extrapolation of the spectra obtained by the 
analysis of LIGO data in the frequency band 40 Hz < / < 500 Hz. The lower dotted line corresponds 
to a 95% upper limit of Ho < 6.9 x 10“ 6 , while the upper dotted line corresponds to an upper limit of 
H 3 < 0.0052 at a reference frequency of 900 Hz. 

on the previously published upper limits. 

We can also compare our results with indirect upper limits on the stochastic gravitational wave 
background. In this band, the most stringent constraints come from Big Bang nucleosynthesis 
(BBN) and measurements of the cosmic microwave background (CMB). The BBN bound con- 
strains the integrated energy density of gravitational waves over frequencies above 10 -10 Hz, based 
on observations of different relative abundances of light nuclei today. The BBN upper limit is [6] 


where N v is the effective number of neutrino species at the time of BBN. The most recent limit on 
this number is ( N v — 3) < 4.6 x 1CT 2 , given by measurements of the relative abundances of light 
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elements and taking into account neutrino oscillations [30]. The CMB limit also constrains the 
integrated gravitational wave energy density, and is obtained from the observed CMB and matter 
power spectra, as these would be altered if there were a higher gravitational wave energy density 
at the time of decoupling. The CMB upper limit [31] is 


Our upper limit is not sensitive enough to improve on these indirect upper limits, however, these 
indirect bounds only apply to a background of cosmological origin, whereas the bound presented 
here applies to astrophysical signals as well. 

VI. CONCLUSIONS 

Data acquired by the LIGO and Virgo interferometers have been analysed to search for a 
stochastic background of gravitational waves. This is the first time that data from LIGO and Virgo 
have been used jointly for such a search, and we have demonstrated that the addition of Virgo 
increases the sensitivity of the search significantly, reducing the error bar by 23% even though 
the length of time for which Virgo was taking data was approximately one fifth of the time of 
the LIGO run. The upper limit obtained with the LIGO interferometers only is the most sensitive 
direct result in this frequency band to date, improving on the previous best limit, set with the joint 
analysis of ALLEGRO and LIGO data, by a factor of ~ 7. 

Adding Virgo improves the sensitivity across the frequency band, largely due to the the addition 
of pairs which have different overlap reduction functions. This enables us to cover the frequency 
band more evenly, as well as effectively increasing the total observation time. We can see that the 
sensitivity of the search is much improved by adding Virgo by comparing the standard deviations 
in Table III. However, in this case, the increased sensitivity did not lead to a decreased upper limit, 
as the joint estimator of obtained by the the full LIGO- Virgo search was higher than the the 
estimator obtained by the LIGO-only analysis. 

As part of this analysis, we have also developed a method of marginalizing over the error on the 
amplitude calibration of several interferometers. The methods used in this paper will be useful for 
future analyses of data from the network of interferometers, which we expect to grow, eventually 
including not only interferometers in North America and Europe, but also hopefully around the 
world. 
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